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Abstract. We begin the construction and the analysis of nonosciUatory shock capturing 
methods for the approximation of hyperbolic conservation laws. These schemes 
share many desirable properties with total variation diminishing schemes, 
but TVD schemes have at most first order accuracy, in the sense of 
truncation error, at extrema of the solution. In this paper we construct 
a uniformly second order approximation, which is nonosdllatory in the sense 
that the number of extrema of the discrete solution is not increasing in time. 

This is achieved via a nonosdllatory piecewise linear reconstruction of the 
solution from its cell averages, time evolution through an approximate 
solution of the resulting initial value problem, and averaging of this 
approximate solution over each cell. 
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1. Introduction. In this paper we consider numerical approximations to weak 
solutions of the scalar initial value problem (IVP) 

(1-la) U( + f(u) x = u, + a(u) u x = 0 

(1-lb) u(x,0) = «o(x) . 

The initial data t*o(x) are assumed to be piecewise - smooth functions that 
are either periodic or of compact support. 

Let vf = v h (xj, t„), xj - jh, t n - nr, denote a numerical approximation in 
conservation form 

(1.2a) v/ + l = v; - (E, ■ v”)j 

Here E h is the numerical solution operator, X = r/h, and f)+\n, the 

numerical flux is a function of 2k variables 

«: 

C 1 - 2 * 5 ) fj+in = f(yf-k*uSj+k) 

which is consistent with (1.1a) in the sense that 

(1.2a) f(u,u,...,u) = /(u) . 

We consider the numerical approximation v h (x,t) in (1.2) to be a 
picccwise-constant function 

(1.3) v A (x,r) = vf, Xj-y 2 < x < x j+y2 , nr < t s: (n + 1 )t . 

Accordingly we define its total variation in x to be 

(1.4) TV(v») = TV(v h (-,t„)) = 2 | vf +1 - vf\ . 

J 
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If the total variation of the numerical solution is uniformly bounded in h for 

Osrsr 

(1-5) 2V(v*(-.f)) * C • rv(uo) 

then any refinement sequence h - 0, r = 0(h) has a subsequence hj - 0 so 
that 

(1.6) v h] — > u 
where u is a weak solution of (1.1). 

If all limit solutions (1.6) of the numerical solution (1.2) satisfy an entropy 
condition that implies uniqueness of the IVP (1.1), then the numerical scheme is 
convergent (see e.g [3], [12]). 

Recently we have introduced the notion of total variation diminishing (TVD) 
schemes (see [3]), where the approximate solution operator is required to dimin- 
ish the total variation (1.4) of the numerical solution at each time-step 

(1.7) 7V(v" +1 ) 2 S TV(y n ) ; 

these schemes trivially satisfy (1.5) with C = 1. Some early work along these 
lines was done by van Leer in [15]. 

TVD schemes are non-osdllatory in the sense that the number of local 
extrema in the numerical solution is diminishing in time (as is customary we use 
" diminishin g" loosely as short for "non-increasing", throughout this paper). 
Moreover, the value of an isolated local maximum may only decrease in time, 
while that of a local minimum may only increase. 

We were able to construct TVD schemes that in the sense of local truncation 
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error are high-order accurate everywhere except at local extrema where they 
necessarily degenerate into first-order accuracy (see [4], [13], [10], [11], [14]). 
The perpetual damping of local extrema determines the cumulative global error of 
the "high-order TVD schemes" to be 0(h) in the L„ norm, 0(1 r^ 2 ) in the L 2 
norm and 0(h?) in the L* 'norm (see [17]). 

hi this paper we introduce a larger class of non-osdllatory schemes, in which 
the solution operator is only required to diminish the number of local extrema in 
the numerical solution. Unlike TVD schemes, which are a subset of this class, 
non-osdllatory schemes are not required to damp the values of each local 
extremum at every single time-step, but are allowed to occasionally accentuate a 
local extremum. 

In a sequence of papers, of which the present paper is the first, we show how 
to construct non-osdllatory schemes that axe uniformly high-order accurate (in 
the sense of global error for smooth solutions of (1.1)). In this first paper we 
describe a second-order accurate scheme of this type. 

The fact that the number of local extrema in the numerical solution may only 
diminish in time is suffident by itself to guarantee that the application of the 
scheme to monotone data results in a monotone function. Thus non-osdllatory 
schemes, like TVD schemes, are monotonitity preserving. In particular, when 
applied to a step-function, they do not generate spurious oscillations. 

We note that since the number of local extrem a in the solution of non- 
osdllatory schemes is bounded by that of the initial data, uniform boundedness of 
its total variation (1.5) follows immediately if the maximum norm of the solution 
is shown to be uniformly bounded. 
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2. Design Principle and Overview 

In this section we describe how to construct a non-osdllatory scheme that is 
uniformly second-order accurate. 

Integrating the partial differential equation (1.1a) over the computational cell 
{Xj-Vb x j+in) x (*n, *«+i) we get 

(2.1a) 5T 1 = i 7 - 1 / 2 OO " h- wMl 

where 

(2* lb) fj+V2(u) = ~ 0) dt > 

and 

( 2 * lc ) “? = \ i x ^ «(*>0 dx . 

iJ 

We observe that although (2.1a) is a relation between the cell-averages uj 
and uj +1 , the evaluation of the fluxes /j+i^(u) in (2.1b) requires knowledge 
of the solution itself and not its cell-averages. 

As in Godunov’s scheme and its second-order extension by van Leer [16] and 
Colella and Woodward [2], we derive our scheme as a direct approximation to 

(2.1) . We denote by vj the numerical approximation to the cell-averages uj of 
the exact solution in (2.1c), and set vj 5 to be the cell-averages of the initial data. 
Given v" = {vj 1 } we compute v fl+1 as follows: 

First we reconstruct u(x,t n ) out of its approximate cell-averages {vj*} to the 
appropriate accuracy and denote the result by L(x; v"). Next we solve the IVP 

(2.2) v, + f(v) x = 0, v(x,0) = L(x; v n ) 
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and denote its solution by v(x,f). Finally we obtain vj +1 by taking cell- 
averages of v(x,t) 

(2.3) v* +1 = ± j [ I '*“ v(x,7) ix . 

J h x s-\a 

The averaging operator in (2.3) is non-osdllatory, therefore the number of 
local extrema in v rt+1 (interpreted as a mesh-function or the piecewise-constant 
function (1.3)) does not exceed that of v(x,t). Assuming v(x,f) to be the exact 
solution of (2.2) implies (since the exact solution operator is TVD) that the 
number of local extrema in v(x,t) is less than or equal to that of 
v(x,0) = L(x; v"). Therefore if the number of local extrema in L(x; v") does 
not exceed that of v n , then the resulting scheme is non-osdllatory. 

We condude that the design of non-osdllatory high order accurate schemes 
essentially boils down to a problem on thejevel of approximation of functions: 
Given cell-averages uj of a piecewise- smooth function u(x), reconstruct u(x) 
to a desired accuracy. Prior to studying this problem we tackle another related 
question in approximation of functions, that' of constructing a non-osdllatory 
high-order accurate interpolation of piecewise-smooth functions. 

In section 3 we construct a non-osdllatory piecewise-parabolic function 
Q(x; u ) that interpolates a piecewise-smooth function u(x) at the mesh points 

(2.4a) Q(xj; u) = u(xj) 

and satisfies, wherever u(x) is smooth, 

(2-4b) Q(x;u) = u(x) + 0(h 3 ) 

( 2 -4c) -£ Q(x ± 0; u) = w(x) + Oih 2 ) . 
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In section 4 we make use of this non-osdllatory piecewise- parabolic interpo- 
lant to design a non-osdllatory reconstruction of a piece wise-smooth function 
from its cell-averages. As in [16], [2], [5], and [9] we take L(x; u) to be the fol- 
lowing piece wise-linear function 

(2.5a) L(x; u) = uj + Sj(x - Xj)/h for [x - x ; -| < hi 2 , 


Unlike the above references that present "second-order accurate" TVD 
schemes, we compute the slopes SJh from Q(x; u) by 


(2.5b) 


Sj/h = m 


-£Q(xj-0;u), 


dx 


0 \u) 


Here m(x,y) is the min mod function 

(s • min(Jx|,ty|) if sgn(x) = sgn(y) = j 
(2.6) *n(x,y) |q otherwise 

We show in section 4 that L(x; u) is a proper reconstruction of u(x) in the 
sense that whenever u(x) is smooth 

(2.7a) L(x; u) = u(x) + O^ 2 ) 

and 

(2.7b) L(x; S) = u(x) + 0(h 3 ) , 

Here u(x) = A -1 f 1 ^ 2 u(x y) dy and L(x; u) — 

/t/2 

= A -1 f L(x + y;u) dy; like Q(x; u ), the latter is also a non-osdllatory 

h/2 

piecewise- parabolic interpolant of u(x), 


(2.7c) 


L(xj; u) = u(xj ) . 



t 


We remark that the "second-order accurate" TVD schemes described in the 
above mentioned references use a slope Sj/h in (2.5a) that approximates 
(d/dx) u(xj) to 0(h), and their loss of second-order accuracy at local extrema 
points is due to lack of smoothness of the coefficient in the 0(h) term at these 
points 1 . This problem is circumvented in the present scheme by talcin g s/h to 
be (2.5b) which is an 0(h?) approximation to (d/dx) u(xj). Unfortunately there 
is a price to pay for this extra accuracy, namely the loss of the TVD property. As 
in TVD schemes 

(2.8) 7V(v n+1 ) as 7V(L(-; v*)) , 

however here 

TV(L('\ v")) a TV(v n ) 

and indeed the scheme may occasionally increase the variation of the numerical 
solution. Although we prove that the scheme is non-osdllatory we have not been 
able as yet to complete a proof of unif orm boundedness of the total variation of 
the numerical solution; this is due to lack of techniques to verify uniform bound- 
edness of the maximum norm of the numerical solution. 

In section 5 we study the proposed scheme in the constant coefficient case. 
We verify that it is uniformly second-order accurate, examine its behavior at local 
extrema points and get estimates for the possible increase in total variation per 
time-step. 

In this paper where we consider numerical schemes of the form (1.2) that are 

1. We repeat that the results of [8] and [11] imply that TVD schemes, no matter how 
they are constructed, must have this loss of accuracy at local e xtr e m a 
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derived from approximating the relation (2.1), it is only natural to consider trun- 
cation error in the sense of cell-averages, i.e. we say that the scheme (1.2) is 
second-order accurate if 

(2.9) S” +1 = E h • a" + 0{h 3 ) 

where u" is the cell-average (2.1c) of the exact solution. Since 

(2.10) u(x) = u(x) + 0{h?) 

whenever u(x) is smooth, (2.9) holds also for pointwise values of the solution. 
However, in the context of 3rd and higher order accurate schemes, this difference 
in definitions of truncation error will be not only conceptual but of practical 
importance as well. 

Up to this point we have assumed that v(x,r) in (2.3) is the exact solution 
to (2.2). The resulting scheme 

(2- Ha) v/ +1 = vj - X[/ y+1/2 (v) - /;_i^(v)] , 

where fj+\n(y) is (2.1b) applied to v(x,t), 

(2- lib) fj+in( v ) = J JJ T f(y(x,t)) dt , 

is certainly second-order accurate in the sense of (2.9). Starting with the exact 
cell-averages vf =» uj in (2.11) we get from (2.7a) that 

(2.12a) v(x,t) — u(x,t + t n ) + 0(h*) for Osrsr 

and consequently 

(2.12b) j) + V2(v) = / J + V2(u) + 0(*2) , 

which implies (2.9) due to the sufficient smoothness of the coefficient in the 
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0(hr) term in (2. 12b) . 

In section 6 we replace the exact solution v(x,r) in (2.3) by an approximate 
one, which we denote by v„(x,t) , This approximate solution is conservative. 
TVD, and second-order accurate in the sense of (2.12a). Thus replacing v(x,f) 
in (2.3) by this approximate solution results in a conservative scheme that is non- 
osdllatory and uniformly second order accurate. 

We remark that an alternative approach to the above is to approximate 
fj+uiiy) is (2.11b) by using a midpoint rule (or trapezoidal rule) for the 
integral and by replacing v(x,r) with a non-oscillating second-order accurate 
approximate one v„(x,r) (see [16] and [2]). The resulting scheme 

(2.13a) v/ +1 = v/ - - fj.vd 

(2.13b) fj+ m = /(v„(* /+I/2 , r/2)) 

is certainly second-order accurate, and it isnon-osdllatory in the constant coeffi- 
cient case. Since we have not used the cell-averaging (2.3) to derive this scheme, 
we cannot ascertain in general that the resulting scheme is non-osdllatory. 
Nevertheless, our numerical experiments as well as many other experiments in 
the context of TVD schemes (see e.g. [1], [2]) demonstrate that the numerical 
results are non-osdllatory in many (if not all) applications. 

In section 7 we present some numerical experiments that compare the present 
scheme with a typical "second-order accurate" TVD scheme. 

3. Nonostillatory interpolation. 

The oscillatory nature of second order accurate Lax-Wendroff type schemes 
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results from a Gibbs phenomenon associated with hi gh -order interpolation across 
discontinuities, hi this section, as a preparatory step towards designing a nonos- 
cillatory approximation to (1.1.), we construct a non-osdllatory piecewise- 
parabolic interpoiant Q(x; u) to a piecewise-smooth function u(x) such that 

( 3 * la) Q(x { ; u) = u(x,) 

(3.1b) Q(x; u) = q i+ y 2 (x; u ), Xi s; x 2 S x i+1 , 

where $ i+1/2 is a quadratic polynomial, and 

(3.1c) Q(x; u) - u(x) = 0(h\ Q(x ± 0; u) - u(x) = 0(h *) 
wherever u(x) is smooth. 

Q(x; u) is non-osdllatory in the sense that the number of its local extrema 
does not exceed that of u(x). 

Since 

4 i+ 1/2C*;; u ) = u h ft+iafo+i; “) = “f+i 

it can be written in the form 

(3.2a )q i+ y 2 (x-ji) = «,• + d i+y2 u '( x “ x i )/* + y°»+l/2 u (x - x,)(x - x i+1 )/h 2 

where 

(3.2b) di+v2 “ = «i+i “ “« 

and D i+V2 u is yet to be determined. 

Dt+V2 u = h 2 ql’+ i^(x; u), x, +1 . 
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We consider as candidates for 4, -1.1/2 the two quadratic polynomials q { 
and q i+1 , interpolating u( x) at (x,-!, x it x 1+1 ) and (x f , x i+i , x,+ 2), respec- 
tively, and choose qi+y 2 to be the one that is least oscillatory in [x„ x l+1 ]. 
Both qj, j = i and j = i + 1, can be written as (3.2a) with D iJryi u = Dm 
where 


(3.2c) DjU = dj+yjM ~ dj- y 2 u = U j+1 - 2 Uj + Uy_-, . 

Since the least oscillatory of q t and ^, + i can be characterized as the one that 
deviates the least from the line connecting (x,,m,) with (x, +1 ,u J+1 ) we choose 
D i+yl “ in (3.2a) to be 

(3.2d) D i+y 2 u = m(Z) J «, 2?, +1 u) , 

where m(x,y) is the min mod function 


(3.3) 


m(x,y) 


s • min(|x|,bl) ,if sgn(x) = sgn(y) = j 
0 'otherwise. 


If u(x) is smooth in [xj-ijj+i], then qj as a quadratic interpolant of u 
satisfies 

(3-4) ^(x) - u(x) = 0(h 2 ), qj(x ) - u(x) = 0(h 2 ), x j . 1 S x 2S x j + 1 . 

If D,u • Dj+iu a 0 then ^,+1/2 is either ^ or q^y Otherwise we set 
D i+ ^ u = 0, but then smoothness of u implies that Dju — 0(h 2 ) and conse- 
quently ?,+i/2 - qj = 0(A 3 ) for 7 = 1, / + 1. Thus (3.1c) follows from (3.4). 

We turn now to prove that Q(x; u) is a nonosdllatory interpolant of u, 
i.e. that the number of its local extrema does not exceed that of u. We do so by 
showing a one-to-one correspondence between local extrema of Q to those of the 
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mesh function {uj}, the number of which certainly does not exceed that of u(x). 

Q may have a local extremum in either the interior of some interval 
(x,,x i+1 ) or at a mesh point x,. The first case, which will be refered to as 
interior-extremum, occurs when there is a point x* , x. : < x* < x i+i , such that 


-£■ Q(x*‘, ») = 0, but — <?i+ly2 ^ 0 . 

From (3.2a) if follows that Q has an interior-extremum, in (x„ x, +1 ) if and 
only if 

(3.5) l^i+i/2 “I > 2|</ f+1/ 2 u\ . 


ql+ 1/2 = <?,+ i/2CO» the value of the interior-extremum is then 


(3.6) 


«f+l/2 = 


J °/+l/2“ 


^i+1/2 u 



2 ’ 




if D i+ y 2 u < 0 it is a local maximum; if Z> i+1/ 2 u > 0 it is a local minimum. 
Since D i+1/2 “ — »»(D,m, D i+1 u), (3.5) holds if and only if 

(3.7a) Di u ■ D i+ \u > 0 


(3.7b) | Dj uj > 2|d,- +1 /2 u\, j = /,/ + !. 

This implies that ^/+i /2 has a local extremum in (x,, x 1+1 ) if and only if both q, 
and also have a local extremum in (x,, x, +1 ) and of the same kind. Since 
a parabola has at most one local extremum, it follows then that q t does not have 
a local extremum in (x,_i, x,) and q i+ i does not have one in (x l+1 , x, +2 ). 
Consequently Q is monotone in both (xf_i, x,) and (x l+1 , x 1+2 ), but in an 
opposite sense, i.e. u ‘ d i+ u < 0; the latter implies that u has a local 

extremum in [x,-, x J+1 ] and that either u { or u, +1 is a local extremum of the 
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is counted as a 


mesh function {u ; } (for obvious reasons the case u,- = » i+1 
single-extremum). The above analysis also shows that interior-extrema are iso- 
lated, i.e. if Q has an interior-extremum in (x i} x, +1 ), then it is the only local 
extremum of Q in (x,-_ lf x, +2 ). 

We turn now to examine the case that Q has a local extremum at a mesh 
point x,; this will be refered to as a mesh-extremum. The above observation 
that interior extrema are isolated excludes the possibility that Q has an interior- 
extremum in either (xi-i, x,) or (x,, x i+1 ) and consequently Q is monotone in 
these intervals. This implies that “ ' d i+i /2 u < 0 and therefore u,- is a 
local extremum of the mesh function {uj}. This concludes the proof that Q(x; u) 
is a non-osdllatory interpolant of u. 

We next express the non-osdllatory nature of Q in terms of total variation. 
If \Dj+ 1/2 u\ s 2|<fy + 1/2 u| then (2.5) implies that Q is monotone in [x^, x ; >j. 
Thus 

(3.8a) p j+ j /2 m| 2 jd ;+ i /2 “i =► TV[ Xj ^(Q) = ty+i /2 »l - 

If 1 Dj+yz “I > 2j dj+y 2 »j then Q has a local extremum in (x,, x J+1 ) and 
TV [X) , XJ ^(Q) = kj+V2 ~ “/I + l u ;+l ~ Qj+vA • 

Using (2.6) we get 

(3.8b) \Dj+y2 “I > 2|d / +1/2 u| => TV[ Xj Xj ^(Q) 

= \d i+V2 «| + ]p J+i/2 «l [ 3^7 - }] • 

We condude that 


- 15 - 



(3.8c) 0 =s TV(Q) - 2 14+1/2 »| s 2 |0»+V2 «i f TT^ 

j mtM l °m + V2 “ 

s ■7 2 lAn+1/2 “i • 

The sum in the RHS of (3.8c) is taken over the set of indices M of intervals 
(x m , a^ +1 ) in which p m+y2 “I > l¥m+\n “l» i.e. w h crc Q has interior- 
extrema. 

Next we show that if u(x) is a piecewise- smooth function of bounded varia- 
tion, then 

(3.9) lim 7V(Q(-; «)) = TV(u) . 

A-nJ 

We observe that in this case the number of intervals in M is finite and is uni- 
formly bounded by the number of local extrema in u(x). Hence (3.9) will follow 
if we prove that D m +y 2 - 0 as h - 0 for all m € M. To accomplish this we 
show that for h sufficiently small M does not include intervals (*„, x^+j) in 
which u(x) is discontinuous. To see that let us examine the case where u(x) 
has a discontinuity at x € (x„ x <+ i) . Clearly 2 u approaches the size of 
the jump in u while di-y^u approaches zero as h - 0. Consequently 

(3.10a) |D,u/d<j.i /2 **| — |l — d,— y2 u /di+y2 M l ■* 1 as h — 0 

Hence for h sufficiently small 

(3.10b) 2| dt+n u| > pi u\ 2: pi+m A 

which implies i g M. 
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4. Non -oscillatory Reconstruction. 

Let u(x) be a piecewise- smooth function and denote by u(x) its mean over 
(x - h/2jc + h/2), i.e. 

C 4 - 1 ) “(*) “ J «(y) <*y • 

We denote uj — u(x ; ) and refer to these values as cell-averages of u(x). Given 
{m}}, the task in hand is to reconstruct u(x) to 0(h *) in a non-osdllatory way; 
denote the approximately reconstructed function by L(x; u). To achieve 0(h 2 ) 
accuracy it is sufficient to consider L(x; u) to be a piecewise-linear function. To 
make L(x; u) a non-osdllatory approximation we use the non-osdllatory piece- 
wise parabolic interpolation Q(x; u) to compute its slopes as follows: 

(4.2a) L(x; u) - uj + Sj(x - Xj)/h for [x - xj[ < ^ h 

i* 

(4.2b) Sj = h • Q(xj - 0; 5), ■— Q( Xj + 0; 5) j . 

Here m is the min mod function (3.3); <fj+ v 2 “ and D j+1/2 u are (3.2b) and 
(3.2d), respectively. 

We note that L(x; u) may be discontinuous at {xj +1/2 } and that 
(4.3a) L(xf, u) = uj . 

To see that wherever u(x) is smooth 
(4.3b) L(x,u) - u(x) = O^ 2 ) 

we observe that in this case 

(4.4a) u(x) = u(x ) + ©(A 2 ) 
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and therefore it follows from (3.1c) that 

(4.4b) J Sj = u(zj) + 0(A 2 ) = u(x,) + 0(A J ) . 

Consequently the RHS of (4.2a) can be expanded as 
(4.4c) L(x; u) = uj + (x - xj) ~ u(x J ) + C^ft 2 ) 

= u(x) + O^ 2 ) for Jx ~ *j\ < \ h 
and thus (4.3b) follows. 

Denote by L(x; u) the mean value of L(x; u) in (x - ft/2, x + ft/2), i.e. 

(4.5a) £(*; 5) - } £ + “ L(y; a) <fy . 

Using (4.2a) to evaluate the integral in (4.5a) we find 
(4.5b) L(x; u) = + d j+m u 

• (x - xj)/h + (V2)(Sj+i - Sj)(x - xj)(x - x j+1 )/h 2 , 

for xj s: x as x ;+1 

(4-5c) L(x ; ; u) = uj . 

Hence L(x; u), like Q(x; a), is a piecewise-parabolic interpolant of u(x). 
Comparing (4.5b) with (3.2) we find that for xj ^ x ^ x 7 - +1 

(4.6a) L(x; u) - Q(x; u) = y(5 ;+1 - 5/ - D j+y2 »)(* ~ *j)(* “ x ;+ i)/ft 2 . 
From (4.4b) we see that * - * £ + OC* 3 ) (Note that this is true 
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also at local extrema points) and therefore 

Sj + 1 - Sj = h 1 —■ u(xj+y2) + 0(h 3 ) . 

On the other hand (3.2) shows that 

°j+i/2 » = * —2 m) + 0(h 3 ) . 

dx *■ 

Therefore 

(4.6b) Sj + 1 — Sj — Dj+i /2 u — 0(h 3 ) 

which shows that RHS of (4.6a) is 0(h 3 ). Since (2.1c) shows that 

Q(x: u) - u(x) - 0(h 3 ) 
we conclude from (4,6a)-(4.6b) that 
(4.6c) L(x; u) - u(x) = 0(h 3 ) . 

We turn now to prove that L(x\ iT) is a non-osdllatory approximation to 
u(x); this certainly implies that L(x;u) is a non-osdllatory approximation to 
“(*)• We shall do so by showing that TV[ X]tXj ^(L(-;u)), the total-variation of 
L(x\ u) in [xj,xj +1 ], which has the value 

(4.7a) TV^jdOja)) =. i(|S,| + fy+J) + | d j+u2 u - |(S y + S /+1 )| , 
can also be expressed as 

(4.7b) TV [xj ^(L(-; 5)) = max(j d J+1/2 *1) - 

Then it follows immediately form (4.7b), (4.3a) and (3.8) that I is monotone in 
[xj, Xj+i] if and only if Q is; consequently L is a non-osdllatory approximation 
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to u(x) in exactly the same sense as Q is to die interpolated function (see sec- 
tion 3). 

Next let us denote 

(4.8a) - Sf = A • Q(xj ± 0; 5) , 


i.e. 

(4.8b) Sj = d/_i /2 u + y Dj-y 2 u, S* = dj +y2 u ~ y &j+vi u » 
and observe that (4.2b) implies that 

(4.8c) |([S;i + ft +1 |) = |[|m(5-, J/)| + KSf +1 , S/ +1 )|] 

^ ! + i ^/'+ ll ) = Y [ l ^/+]/2 U ~ ~2 ® j + U . 2 “1 + M / + 1/2 M + Y ®/+ l /2 “1 


- max nd,- 4 . ^ u], yl^+ 1/2 “I 


We note that if |d^ +1 /2 “1 ^ l/2jf>/+i/2 **1 ^ en 

• sgn^+^S) 

; 

which in turn implies 


1 


0 


sgn(5 ; ) - sgn(dj+i /2 u) a 0 , sgn(^ +1 ) * sgn(d , +1/2 5) a 0 

It follows then from (4.8c) that the RHS of (4.7a) is \dj + 1/2 “I* This shows that 

(4.9a) \d ]+1/2 q a i |D , +1/2 5| =* 2 V m ^j(L(-; u)) = |</, +u2 q . 
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To complete the verification of (4.7b) we still have to show that 
(4.9b) | dj+n a] < { | D ]+m 21 = 7V Mjrt] (t(-i 5)) = | ID,.! 21 . 


First we observe that 

(4.10) Sf ~ S~ = (</,-+ 1/2 » ~ y f>i+i/2 «) “ (4i- 1/2 “ + Y -Of— 1/2 “) 


- 1 - 

= Di U - —(Dj+ia £ + Oj-l/2 **) 

Since (3.2d) implies 

(4.11) |d,21*A(1d,_ w «1 + 1 d (+w 2D 
we conclude from (4.10) that 

(4.12) (ST - Sf) ' sgn(D, u) * 0 . 

We turn now to prove (4.9b). First let us consider the case that Q(x; u ) 
has a local maximum in (xj, x ;+1 ), i.e. Dj u < 0, Dy +1 u < 0, and 

\ d j+v 2 “1 < y )Pj+ 1/2 “ 1 * 

It follows from (4.12) that 

(4,13a) Sf ^ Sf = dj+i /2 u ~ ~2 ty+i/2 u > 0 

(4.13b) 0 > dj+ 1/2 u + y Z> y -+i/2 m = S/+i 2: S/+i • 

The relations (4.13) and the definitions (4.8a), (4.2b) imply that 
(4.14a) Sj = Sf = d j+y2 « “ Y D j+v 2 " 
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(4.14b) 


Sj+i ~ s j + 1 “ d j+v 2 “ + ~2 D j+yi “ • 

The same analysis shows that (4,14) holds also for the case that Q(x; u) has a 
local minimum in (xj, Xj+{). (4,9b) follows immediately from (4.14) and (4.7a). 

We note that since L(x;u) is continuous at xj 
(4.15) TV<M-. 5)) = 2 TV[ S)) = 2 ma*L +u2 «|, i |D J+ m S| 

J i " 

V / 

= 2 |dj- + i/2 m] + 2 4- Pm+\J2 - Ki + 1/2 «*l) • 

j m£M (4 J 

Here M is the set of indices of intervals (x mt x m ^. i ) in the interior of which L 
(and also Q) has a local extremum. The number of these intervals is finite and 
is bounded by the number of local extrema of u(x). Comparing (4.9) with (3.8) 
we note that 

<« 

(4-16) 7V(L(;5))a2V(G(;5)). 

5. The constant coefficient case. 

In this section we study the constant coefficient case 

(5.1) u, + au x = 0, a = const . 

The exact solution of the IVP (2.2) is 

(5.2) v(x,t) = L(x - at; v n ) . 

Hence our scheme (2.3) is 

(5.3a) vj +1 = — 1+Vi L(x - ar; v fl ) dx = L(xj - ar; v n ) 
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where L is (4.5a). We have shown in section 3 that the number of local 
extrema in L(x ; v'*) does not exceed that of v n . Since v"'*' 1 in (5.3a) is a cell- 
average of L, it follows that the number of local extrema in v n+1 does not 
exceed that of v n , and consequently the scheme (5.3a) is non-osdllatory. 

Using (4.5b) in (5.3a) we get the following expression for the scheme 

(5.3b) v? +1 = (E„ ■ v")j 

' 

Vf ~ My- V 2 v" - 1/2 Ji(l - iL)(SJ ~ 5;_!) if a > 0 

= v; - M/+ 1/2 v" + 1/2 |i(l + m.)(5; +1 - Sf) if a < 0 ; 

£/, denotes the operator form of the finite difference scheme; \l - \a, the 

CFL-number, is assumed to satisfy 

(5.3c) IM'I • 

We turn now to prove that (5.3) is second order accurate in the sense of 
(2.9), i.e. if u"(x) denotes the mean value (4.1) of u(x,t„) then 

(5.4a) SJ* 1 - (E„ ■ i 7)j = 0(A 3 ) . 

To show that we observe that m the constant coefficient case (5.1) 

= iP(xj ~ ar), and by (5.3a) (E h * iT)j = L(xj - at; u"). Hence the LHS 
of (5.4a) is nothing but 

(5.4b) iT(xj — ar) — L(xj — ar",u n ) , 

which is 0(h 3 ) as a direct consequence of (4.bc). 

Next we study the time-dependence of the total variation and the maximum 
norm of the numerical solution (5.3). In section 2 we have pointed out that 
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(5.5a) 7V(v' ,+1 ) s TV{L{ \ v")) . 

Using (4.15) and (5.5a) we get the following upper bound on the possible growth 
of the total variation of the numerical solution per time-step 

(5.5b) JV(v” +1 ) - 7Y(v") s 2 ■ 

Here Af„ is the set of indices of intervals (x m , x m+1 ) in the interior of 
which L(x; v fl ) (and also Q(x; v")) has a local extremum. The number of these 
intervals is finite and remains uniformly bounded in time by the number of local 
extrema in the initial data. 

Gearly the upper bound (5.5b) is overly pessimistic. It estimates the possi- 
ble increase in variation in the reconstruction step due to replacing the cell- 
averages vf by the piecewise-linear function L(x; v n ). It does not take into 
account the possible decrease in variation in the averaging step (2.3), resulting 
from doing just the opposite, i.e. replacing the piecewise-linear function 
L(x - ar\ v 71 ) in (5.2) by its cell-averages (5.3a). 

In the following we shall examine the temporal behaviour of the local 
extrema of the numerical solution and its total variation by analysing the explicit 
values of the cell-averages vj* +1 given by (5.3b). To simplify our presentation 
let us assume that a > 0. 

First we note that (4.8b) implies 

(5.6a) ]Sj - Sj-d = |Sy| + ft., | a 2 maxfl, v*|, |D;- W v"l) . 

Hence 
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(5.6b) Uj-vi v”! a j Pj-V2 v"! =» hf>-wl = IS? - v*| s 2 . 

Rewriting (5.3) in this case as 

(5.7a) vj* +1 = vf - a dj v n - V2 y.(l - *jl)7 ; _i/ 2 dy-i/2 y " 


= (1 - (Tj-i/2) vf + yf-i 

where 


(5.7b) 07 - 1/2 = M- + J M*(l ” pO*fy-l /2 • 

We see that the CFL condition 0 < p. :s 1 and (5.6b) imply that 
(5.7c) 0 ^ (Tj-yi 55 1 5 

thus we conclude 


(5.8) 


to-1/2 v”| 2 i |D ; - w v«-| => V?*' ( [vf.,, vf] 


Relation (5.8) shows that if v" is monotone for Ji ^ j ^ Jr, i.e. 
vj L ^ v Jl+i s • • • s v J k , or vj L a v /i+1 4 ' • a vj K , then v" +1 is mono- 
tone for / £ + 1 s ; s Jr, and in the same sense. Relation (5.8) also shows 
that mesh-extrema of v", i.e. those for which Q has its local e x t rem um at a 
mesh point, are being damped at the n-th time-step. Namely, 


(5.9a) \d J±1/2 t"|2 y Pj ± u2 **l. yf-i 5 v ” 2 => vf/, 1 ) s vf 

(5.9b) ty tW v"| a | |D /±lfl v*|. vf., a vf s vf., => mm(vf +1 , vf?, 1 ) a= vf 


We turn now to consider interior local extrema of v", i.e. those for which 
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Q has its local extremum in the interior of some (*,, x, +1 ). We recall that such 
an extremum is characterized by v fl j < 1/2 \D i+ y 2 v fl j and that Sf and 

5f +1 in this case are given by (4.14); therefore - Sf = Df+i /2 v". From 
(5.3a) and (4.6a) we see that in general 

(5.10a) vf+i 1 - Q(x i+ i - ar; vf = ~ ix(l - p,)(D 1 + y2 v fl - Sf +1 4- Sf) . 

Hence 

(5.10b) | 4 +i/2 v n \ < j 1 D i+ v 2 V| => vf+j 1 = Q(x i+1 - at; v n ) . 

Relation (5.10b) confirms the second order accuracy of the scheme at local 
extrema. Although it does not necessitate accentuation of the extremal values, as 
vf+i 1 in (5.10b) may still be in [vf, vf +1 ], it does allow vf+i 1 to deviate from 
this interval by as much as 

(5.10c) | |D i+ M »*|(H +M v*HDi+m »”| - VI) 1 

Thus (5.10b) is the essential difference between the present scheme and the 
"second order" TVD schemes. 

A similar analysis, which we do not present here, shows that if vf is a 
mesh-extremum then vf +1 , j — i, i + 1, relates to Q(xj - ax; vf in the fol- 
lowing way: 

(5.11a) vf +1 s Q(xj - ar; vf, j = i, i + 1, if vf is a maximum 
(5.11b) vf +1 :S Q(xj - ar; vf, j = i, i + 1, if vf is a minimum . 

From (5.9)-(5.11) we deduce the following relation between the total varia- 
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tion of the numerical solution and that of its piecewise-parabolic interpolant Q: 
(5.12) TV({Q(xj ~ ar; v n )}) s rV(v fl+1 ) ^ TV(Q('\ v*)) . 

The LHS of (5.12) is the totai variation of the mesh function 
{Q(xj - at; v")}. Relation (5.12) suggests to consider an equivalent definition 
TV of the total variation of the numerical approximation of the form 

2Y(v fl ) = 7V(fl(*; v-)) 

with the hope that the scheme (5.3) is TVD with respect to this modified defini- 
tion. Unfortunately our numerical experiments have shown that there are 
instances, although rather rare, that TV(v n ) is increasing with n; the same is 
true for TV(y n ) = IV(I(*; v fl )). 

As we have mentioned in the introduction, because of the nonosdllatory 
nature of the scheme, unif orm total variation boundedness of the numerical solu- 
tion is implied by uniform boundedness of its maximum norm. If we follow a 
particular local maximum of the initial data we see from (5.9)-(5.10a) that it 
actually decreases most of the time, and whenever it does increase (5.10c) and 
(3.10) suggest that it does so by a "small amount" that vanishes with h - 0. 

Since the initial data is only piecewise-smooth we have not been able as yet to 
rigorize these arguments. 

We remark that our numerical experiments clearly indicate that in a normal 
computational situation the maximum norm of the numerical solution is indeed 
uniformly bounded. We feel that our inability to prove this fact stems only from 
lack of theoretical tools to analyse pointwise regularity of the numerical solution. 
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6. The nonlinear case. In this section we describe an approximate solution 
v n (x,t) of [5] for the IVP (2.2) 

(6.1) v t + f(v) x = 0, v(x,0) = L(x; v n ) . 

This approximate solution is consistent with the conservation form of the equation 
(6.1) in the sense that the cell-averaging (2.3) results in a scheme in conservation 
form i.e. 


(6.2) v/ +1 = ^ f x ^ v„(x,t) dx = v] - \(fj+u2 ~ fj- in) 

where the numerical flux fj+m is consistent with f(u) in the sense of (1.2c). 
Furthermore, the approximate solution operator is TVD 

(6.3) 2V(v n (-; /)) s 7V(v n (-; 0)) = 7V(L(*; v")) for 0 =s rs t 


and thus by the reasoning presented in section 2, the resulting scheme (6.2) is 
non-osdllatory. 

We turn now to outline the derivation of this approximate solution. To sim- 
plify our presentation we ignore entropy considerations and refer the reader to 
future papers for details of appropriate modifications. We assign to the point 
Xj+ 1/2 a characteristic speed that corresponds to the Rankine-Hugoniot speed 
njf+i /2 of the two neighbouring cell-averages vj and vj +l 


f /Of+i) - f(yf) 


(6.4a) 


a )+m - 


v/+i - vf 


a (yf) 


if vf * vf +l 


if vJF = v / +1 


and denote by a(x) the piecewise-linear interpolant of i.e. 
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(6.4b) a(x) = a 7 -_y 2 + 0*/+i/2 ~ a j-\n) ’ (* “ *7-1/2)/* 

for Xj-V 2 ^ * 25 *7-1/2 - 

The approximate solution v„(x,t) is defined by specifying its constancy 
along the approximate characteristics 

(6.5a) x(t) = x Q + a(xo) • f 

i.e. 

(6.5b) v„(x(0, 0 * v n (*o»0) = L(x 0 ; v n ) . 

Using (6.5a) and (6.4b) to express x 0 in terms of x and t 
6.5C) X Q (x,f) = Xj-y 2 + h [x - Xj-y2(t)V[Xj+y2(t) ~ *7- l/zto] 

for Xj- m (t) <x< x j+yi (t) 
we get from (6.5b) that 

(6.6a) + * • . v ») 

for xj- m (t) <x< Xj+vfc) 

where 

(6.6b) * i+ i 72(0 = *i+i/2 + * * ®i+l/2 • 

Taking cell-averages of the approximate solution (6.6) we get the conservation 
form (6.2) 

(6-7a) vf 1 = v/ - X(7 ;+1/2 - fj-yj 
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with the numerical flux 


(6 

where 

(6.7c) Sj = S"/[ 1 + M«/+i/ 2 ~ «jf— 1 / 2 )] • 

Note that (6.7) is identical to (5.30) in the constant coefficient case. 

We turn now to prove that the scheme (6.7) is uniformly second-order accu- 
rate in the sense of (2.9). Starting with the exact cell-averages v" = uj in (6.7) 
this amounts to showing that 

(6.8a) fj+yi = fj+\rM + Oih 2 ) 

with a sufficiently smooth coefficient in the OQt 2 ) term; here fj+m. is the 
numerical flux (6.7b) computed with the exact cell-averages, and 1 / 2 ( 1 *) is 
(2.1b). We shall do so in two steps: first we shall show that 

(6.8b) f j+yi (u) = ^[f(L(x j+ Vi; u")) + f(L(xQ(xj+v 2 , t); 0)1 + OCh 2 ) , 

where x 0 (x ;+1/2 ,t) is (6.5c), and then we shall verify that 

(6-8c) yl y(I(x ;+1/2 ; 5")) + /(I(x 0 (x y+ly2 , t); u"))] = fj+V2 + 0(h 2 ) . 

Special attention will be given to the smoothness of the 0(h 2 ) coefficients. 

To show (6.8b) we start by using the trapezoidal rule to approximate the 
integral in (2.1b); we get 

(6.9a) fj j. 1 / 2 ( 1 *) = ~[f(u(xj+ !/ 2 , *„)) 4* f(u(x j+y2 , t „ + t))] + 0(h 2 ) . 


f(vf) + 1/2 a.j+1/2 (1 - X aj-1/2) ■ Sj if <xj. 

/(vf+i) - 1/2 5} + 1 / 2(1 + X a jJr yjj ■ Sj+i if aj. 
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The smoothness of the OQt 2 ) term follows from that of f(u ) and u(jr,f). Next 
we observe that a(x) in (6.4b) approximates a(u(x,t n )) to OQt 2 ), and there- 
fore we can use the approximate characteristic line (6.5c) to trace 
u(xj+ 1 / 2 , t n + t) to u(x Q (xj + y 2 , t), t„) with OQt 2 ) accuracy; consequently 

(6.9b) f(. u ( x j+V2> * n + t)) = f(u(*o(xj+v2> t), O) + 0(h 2 ) . 

Finally we obtain (6.8b) by approximating u(x,t n ) in (6.9a) and (6.9b) to 
OQt 2 ) by L(x; u") (see (4.4)). The smoothness of the OQt 2 ) term in this 
approximation is due to (4.4c); 

Sj = h- «*(*, t n ) + 0(h 2 ) . 

(We recall that the degeneracy to first order accuracy at local extrema points of 
some "second-order accurate" TVD schemes is due to lack of smoothness there of 
the OQt 2 ) term in (2.7a)). 

i« 

We turn now to verify (6.8c). First let. is consider the case 5) +1 ^ 25 
L(*j+v. 2J S’ 1 ) = uj + jSf, 


L(x 0 (x ;+1 ^, r)’,u n ) = uj + 





1 + \(aj±ia - Sj-y£ 


s? 


} 


and expand the LHS of (6.8c) around u!j. We get 

(6.10a) f(uj) + y a(u^)( 1 - X aj+y£ Sj + -|-fl'(u^)[(l - X 5/- 1 / 2) 2 


+ (Kaj+y^^iSj) 2 + 0(h 2 ) = fj+ 1/2 + Y (2 X 2 a 2 - l) a' ■ (u x ) 2 \ J+y2 + OQt 2 ) . 
Similarly in the case a/+i /2 2 0: 
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— 1 

L(Xj+U2>'* n ) = u j + 1 - J S Ul> L ( x o( x j+l/2’ t);**") 


Xa 


[2 1 + \(aj+y2 ~ ^f+ 1 / 2 ) 


= «j+i “ 

we expand the LHS of (6.8c) around u^ +1 to get 
(6.10b)/(u£+i) - j a(uj +1 X 1 + \aj+y£Sj+i 


7+1/2 




+ | a'(«7 +1 )[(l + \b} + k) 2 + (XSJ+i^ 2 ] ■ (i/+i) 2 + 0(h 3 ) 


= /j|+i /2 + y (2X.V - 1) • a' • (b.,) 2 [, + V 2 + °(* 3 ) • 

We see from (6.10a) and (6.10b) that independently of the sign of aj +lf 2 , the 
O^ 2 ) term in (6.8c) is the same, namely 

Y (2 \ 2 a 2 - 1) • o' ■ (m x ) 2 |/ + i^ • 

This completes the proof that the scheme (6.7) is second-order accurate in the 
sense of (2.9) wherever u{x,t) is smooth, including local extrema and sonic 
if = 0) points. 

Remarks'. (1) The numerical flux (6.7b) can be rewritten as: 

(6.11) f , +u2 = y[/(v/) + / 0 ?+i) - (v; + i - v/) 


+ max(0, a j+1 x) • (1 - Xa)-i^) • 5/ 

- min(0, oj+vd • (1 + Xa J+3r2 ) * S ; +i| . 
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(2) Our proof that the scheme (6.7) is non-osdUatory is based on the 
representation of (6.7) as the cell-average (6.2) of the non-osdllatory approxi- 
mate solution v n (x,t) in (6.6). To ensure that v„(x,t) remains univalued for 
0 =s t s t we have to restrict the time-step t so that for all j 

(6.12a) x j+y2 (j) > *,-1/200 . 

This implies the condition 

(6.12b) t maxj(aj - ty+ 1 / 2 ) ^ h . 

On the other hand, to derive the particular form of the numerical flux (6.7b) we 
have made the assumption 

(6.13a) 1*,+ 1/2OO “ *7+1/2! < h for all j , 

which implies the condition. 

(6.13b) t max ; |a}j.’i^| < h . 

The two time-step restrictions (6.12b) and (6.13b) are characteristic to 
Godunov- type schemes. The practice of reconciling the two conditions by 

t maXj|av+i/2l 2 hi 2 

is completely unnecessary: The scheme (6.7), as the original Godunov scheme, 
rem ains non-osdllatory (or TVD in the case analysed in [10]) for the full CFL- 
restriction (6.13b). The reasoning for this observation is as follows: The approxi- 
mate solution (6.6) under the CFL restriction (6.13b) may fail to be univalued in 
the ;-th cell only if aj-m > 0 and a<+ip_ < 0. In this case the numerical fluxes 
fj ±\/2 as denned by cell-averaging in the neighboring cells j ± 1, remain (6.7b). 
Thus the only thing that needs to be changed in this case is the definition of 
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v n (x,T ) in the y-th cell. 


(3) We observe that once v„(x,r) is defined globally in (6.6) there is no 
intrinsic need to average it on the original mesh. We may average it on different 
intervals and still conclude that the resulting approximation is non-osdllatory and 
conservative. Furthe r mo re , the construction of the interpolant Q, the approxi- 
mation L and the approximate characteristic field a(x) needed to define 

v n (x, t ), does hot depend on the uniformity of the mesh. Therefore the scheme 
(6.7) generalizes immediately to non-uniform moving meshes. Of particular com- 
putational interest are the self-adjusting moving grids of the type described in 
[12], which make it possible to obtain perfectly resolved shocks and contact 
discontinuities. 

(4) We note that since the approximate solution v n (x,f) in (6.6) is conser- 
vative, it is possible to consider an associated random-choice method obtained by 
replacing the cell-averaging in (6.2) by a sampling with respect to a variable that 
is uniformly distributed in the cell, i.e. 

v? +1 = v.(xj + 8/A,t) 

where 0" is uniformly distributed in [-172,1/2]. Following the reasoning of [7] 
it is clear that the resulting random-choice method is non-osdllatory and that its 
limits are weak solutions of (1.1). Although the random-choice approach has 
many attractive computational features, it has been our experience that in many 
application it is possible to accomplish the same computational goals with a self- 
adjusting moving grid. In this case the use of the latter is preferable as it offers 
gain in resolution without a loss in pointwise accuracy that is assodated with sam- 
pling. 
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7. Numerical Illustration. In this section we compare the new unif ormly 
second-order non-oscillatory scheme of this paper (to be refered to as UN02) to 
the typical second-order TVD scheme (to be refered to as TVD2). Both schemes 
can be written in the form (6.7), i.e. 

(7.1a) v/ +1 = vJF - \(f j+1/2 - fj -j/a) , 


(7-lb )fj+\n 


f(yf) + 1/2 a ;+1/2 ( 1 — + ^(a/4-1/2 ~ a j-ui)\ 

if a j+1/2 ^ 0 

f( v j+ 1 ) - 1/2 5/+ 1 / 2(1 + X.5jf+3^)5"+i/[l + ^(5/+3/2 ~ ^jf+ 1 / 2 )] 

if 4/+1/2 25 0 


(7.1c) 5/ = m(S;,SJ-); 

here a ;+ 1/2 is (6.4a) and m(x,y) is the min mod function (3.3). Sf are dif- 
ferent for TVD2 and UN02: 

iJ 

(7.2) TVD2: Sf ="d j±u2 v n 


(7.3) UN02: 5* = dj±y2 j D j±in y n , 

where dj+ 1/2 and Oj^i /2 are defined in (3.2). 

UN02 and TVD2 are both second-order accurate Godunov-type schemes 
that differ only in the reconstruction step (4.2a): 

(7.4a) L(xyt) = uj + Sj(x - xj)/h for [x - xj\ < h/2 , 

where the slopes of the lines are calculated by (7.3) and (7.2), respectively. 
Therefore we start our comparison on the approximation level. 

In Table 1 and Fig. 1 we present approximations to 
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u(x) = simrx, — 1 ^ x ^ 1. We divide [-1,1] into N equal intervals and 
define 

(7.4b) xj = -1 4- j • 2. , 

The symbols in Fig. 1 denotes values of u, = sinirac ; - for N = 10 in (7.4b). In 
Fig. la we show the piecewise-parabolic interpolant Q(x-ji) (see section 3). In 
Fig. lb we show the piecewise-linear approximation L^^xju) which is (7.4a) 
with (7.1c) and (7.3). In Fig. lc we show the piecewise-linear approximation 
L TVD1 (xyi) which is (7.4a) with (7.1c) and (7.2). We make the following obser- 
vations regarding Fig. 1: (i) Q is a better approximation than L 1 ^ 02 ; L UN02 is a 
better approximation than L TVD2 . (ii) 7V(L UN02 ) > TV(u) > 7V(L TVD2 ). In 
table 1 we quantify the first observation; we list the L„-error and the -error of 
these approximations to sum* for a refinement sequence of N = 10,20,40,80 
in (7.4b). Clearly Q is an 0(h 3 ) a pp ro ximation, while L 1 ^ 02 and L TVD2 are 
0(h 2 ). The error in L 1 ^ 02 is about a 1/3 of the error in L TVD2 . 

In Table 2 and Fig. 2 we present solutions of UN02 and TVD2 for the con- 
stant coefficient case 

(7.6) u t + Ujc = 0, u(x,0) = sin -rrx, — 1 ^ x ^ 1 

with periodic boundary conditions, at t - 2 with t Zh = 0.8. Figs 3a and 3b 
show UN02 and TVD2, respectively, with N — 20 in (7.4b). In table 2 we list 
the Lao-error and -error for a refinement sequence with N - 20,40,80,160. 
Clearly UN02 is second-order accurate in both L« and Li, while TVD2 is 
second-order accurate in Lj but only first order accurate in L*. Fig. 3b demon- 
strates that the degeneracy to first order accuracy at local extrema in the TVD 
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scheme adversely affects the accuracy everywhere (Because the scheme is TVD it 
cannot switch abruptly to second-order accuracy as this would introduce osdlla- 
tious; consequently it takes quite a while to recover the second order accuracy). 

Next we approximate the discontinuous function 


(7.7) 


u(x) = 


-xsin(3'irx 2 /2) , — 1< x < -1/3 

jsin(2irx)| , |x| < 1/3 

2x — 1 — l/6sin(3irx) , 1/3 < x < 1 


which we extend to have period 2 outside [-1,11- 

In Fig 3 we present approximations to <J>(x), using N = 20. Fig 3a shows 
Q(x-,u), Fig 3b shows I^^xjii), and Fig 3c shows L TVD2 (x;u). We again 
observe that Q is better approximation then I 1 ^ 02 , while L 1 ^ 02 is a better 
approximation then L TVD2 . 

In Fig 4 we present solutions of UN02 and TVD 2 for the constant coeffi- 
cient problem (7.6), initial data given by (7.2), and periodic boundary conditions. 
We take t = 2 and r/h = 0.8. Figs 4a and 4b show UN02 and TVD2 respec- 
tively with N — 40. Fig 4b shows the damping effect that the TVD scheme 
imposes due to its degeneracy to first order accuracy at local extrema. 

In Fig 5 we solve the same problems, except we impose boundary conditions. 
At x = -1 we impose the given function (7.7) evaluated at — 1 - t. No boun- 
dary conditions are imposed at x = 1. We implement this numerically using 
UN02 and TVD2 except at the boundary points. There we are in general, unable 
to construct non-osdllatory piecewise parabolic interpolants Q(x,u), so we con- 
struct the only possible parabolic interpolant thru x,-,x 1+1 and the point to either 
the left or right which lies in the region. The analogous procedure is carried out 
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at the reconstruction stage. Figs 5a and 5b again show the results at t = 2 with 
t /h = 0.8. 


The possible introduction of oscillations through the boundary conditions 
does not seem to have degraded the performance of either scheme (in fact the 
opposite is observed). Again the TVD2 scheme shows a damping effect. 

In Table 3 and Fig. 6 we present results for Burgers’ equation 

(7.8) u, + uu x = 0, u(x,0) = a + sin ir(x + 0), -Isis 1 

with periodic boundary conditions and r/h (1 + jaj) = 0.5. The solution to (7.7) 
is smooth for t < 1/ir; at t — Utt it develops shocks. In Table 3a we list the 
Loo-error and -error of UN02 and TVD2 at t — 0.15 for a = 8 = 0 in 
(7.7). This table shows the same asymptotic behaviour as Table 2, except that 
because of the large gradients it shows for a smaller h. 

In Figs 6a and 6b we show results of tJNQ2 and TVD2 for (7.8) with a = 2 
and p = 1 at t — 0.35 with N — 20. In this case the solution to (7.8) 
develops a shock moving with speed 2 beginning at time t = l/-ir ~ 0.318. 

In Table 3b and Figs 6c and 6d we repeat the previous calculations for the 
schemes (2.13): 

( 7 - 9a ) v? +1 = v/ - X(/ ; +i/2 - fj-ird 

(7.9b) f j+m = f(v n (x j+yl , r/2)) = f(L(x 0 (x J+V2 , r/2),v n )) 

where x Q (x ;+1/2 , t/2) is (6.5c), i.e. 


(7.9C tfj+y2 = 


/ ✓ 



1 - A (dj+yi + m s „ 

1 + \(dj + y 2 ~ J , 



1 + A (fly + 3/2 + dj+iri)/2 
1 + \(dj+2/2 ~ dj +h ry)/2 7 , 


if 4; +1/2 s 0 
if &]+v . 2 s 0 


- 38 - 



As we have remarked in section 2, vj 1 ' 1 ’ 1 in (7.9a) is not a cell-average of 
v„(x,r), but only an approximation to it. Therefore it is not necessary to take 
dj +y2 in (7.8c) to be (6.4a). We choose &jj.y 2 so that (7.9c) is continuous at 


a,- +1/2 = 0: 

(7.9d )d j+V2 = 


M+i - - A*? + \ S D 


|W + i - i^-i) - <y? + {*/)! 


We denote the schemes (7.9) with 5" defined by (7.1c) and either (7.2) or 
(7.3) by FVD2 and FN02, respectively. We note that (7.9) is identical to (7.1) in 
the constant coefficient case, and consequently FVD2 and FN02 are nonosdlla- 
tory in the constant coefficient case. Figs 6c and 6d show that FN02 and FVD2 
are also non oscillatory in the case (7.8). Furthermore, Table 3b shows that 
FN02 is much more accurate than UN02 (FVD2 is about the same as TVD2). 

In all previous examples we have presented pointwise calculations; namely, 
we have initialized the numerical solution by taking vj 5 to be the value of the ini- 
tial data at xj, and we have considered vj to be an approximation to u(xj,t n ). 
(Surely this is an acceptable practice for second order accurate schemes.) In Tabie 
3c we repeat the calculation for LJN02 in Table 3a, but now in a sense of cell- 
averages and denote it by AN02. Now we initialize UN 02 for (7,8) with 
a = 0 = 0 by cell-averages of the initial data, i.e. 

(7.10a)vj ) = j- [cos(irx ;+V2 ) “ cosOirx,.^] = ' sin(irx ; ) , 

and regard vj to represent cell-averages of u(x,t n ). To obtain a pointwise 
approximation to u (x,t n ) we first compute point values VJ + ^ of its indefinite 

integral u n {xj+y£ = f ^“(ySn) <ty by 

■*0 
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(7.10b) 


vy+ifl = * i »?. 

*“i'o 

Next we obtain a global piecewise-linear approximation v(x,t„) to u(x,t n ) by 
(7.10c) v(W„) = ^ 2(x-r>) 

where 2 is the piecewise -parabolic interpolant of section 3. Finally we get 
(7.i0d) Q(x t y") = | (v; +1/! - V”_^) = v; . 

Thus the only difference between AN02 in Table 3c, and UN02 in Table 3a 
is the initialization (7.10a), which itself differs only slightly from the mesh values 
of the initial data (since sin(irA/2)/(irA/2) = 1 — 1/6(11 A/2) 2 + 0(A 4 )). 

We remark that cell-averages do play a si gnifican t role when the initial data 

is discontinuous (since they provide information about the location of the discon- 

, ♦ 

tinuity) and in higher-order Godunov-type schemes; this will be described else- 
where. 

We turn now to present calculations with a formal extension of UN02 and 
TVD2 for systems of conservation laws. We consider a Riemann problem for the 
Euler equations of gas dynamics 

(7.11a) u, + f(u) x = 0, w(x,0) = 

(7.11b) u = (p,m,£) r , /(u) = ( m,m 2 /p + P, m(E + P)/ p) 7 " , 

(7.11c) P = (y - 1)(£ - y .vt 2 /p) . 

Here p, m, E and P are the density, momentum, total energy and pressure, 


’u L x < 0 

Ud X > 0 ’ 

V 
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respectively ; we take y = 1.4. 

In the following we apply the extension technique of [3] to UN02 and TVD2. 
The idea is to extend UN02 and TVD2 to systems in such a way that will be 
identical to (7.1) in the scalar case, and will decouple into (5.3) for each of the 
characteristic variables in the constant coefficient system case. To accomplish 
that we use Roe’s averaging for (7.11) (see [13]) 

(7.12a) Y/+ 1/2 = V(yJ, vf +1 ) 

for which 

(7.12b) /(v/ +1 ) - f(yf) = A(v j+1/2 )(v/ +1 - vf), A(u) = df/du , 

and define local characteristic variables with respect to the right-eigenvector sys- 
tem {/?*+i/ 2 }i»i A(vy+i/ 2 ). We extend (6.11) to systems as follows: 

(7.13a) vf +1 = vf - \(fi+v 2 - f } -yd 

(7.13b) f J+y2 = T /(vj 1 ) + /(v/ +1 ) - 2 cf+yz Rf+M 

(7.13 C)c}+V 2 - \<$+V 2 \ d j+V 2 w - ~ ^aj-V 2 )Sj 

+ min(0,a^ +1 ^)(l + 'kaj +y2 )^j+i • 

Here at +iy2 is the 4-th eigenvalue of A(v ;+1/2 ) corresponding to Rj+y2> and 
4+ yjw denotes the component of dj +y2 v = vj +1 - v" in the 4-th characteris- 
tic field, i.e. 

(7.13d) dj+y 2 V = 2 0*/+ V2 w ) R j+ 1/2 • 
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Likewise S* denotes the component of the vector of slopes in the k-th charac- 
teristic field, and is defined as follows: 

(7.13c) sf = + \(aj+i /2 - i 

m(x,y) is the min mod function (3.3). S k ±J are different for VD2 and N02: 

(7.14) TVD2: = d$ ±y2 w 

(7.15) UN02: Sij = d* +vz w 3= | B* ±w w; 

Df± VI" = '"(4+3/2"' “ 4+l/2 w > 4+1/2 w “ 4- 1/2*0 • 

In Figs 7.8 and 7.9 we show numerical solutions of UN02 and TVD2, respec- 
tively, for the Riemann problem (7.1b) with 

U L = (l,0,2.5) r , U R = (0.125,0,0.25) . 

These figures demonstrate that the formal' extension to systems is nonosdllatory 
in this case. Since the solution to the Riemann problem is just constant states 
seperated by waves we do not get to see here the extra resolution power of UNG2, 
except that its numerical solution is somewhat "crisper” than that of TVD2. In 
this calculation we have not employed any artificial compression in the linearly 
degenerate field and therefore the contact discontinuity smears like n 173 , as 
expected. The interested reader is referred to [4], [5] and [10] for a detailed 
description of such compr e ssion techniques, as well as for details of entropy 
enforcement mechanisms. 
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Table 1 : Approximations to u(x) = sin(-irr), with periodic boundary conditions. 


L.-ERROR 


^UN02 


Lr ERROR 




1.420 X 10" 1 


3.558 X 10" 2 


9.163 x 10" 3 



1.494 x 10" 2 i 2.467 x 10" 2 


3.104 x 10” J 7.710 x 10 -4 



9.787 x 10- 


Table 2: Numerical solutions of u, + u x = 0, u(x,0) = sin ttx, -Isis 1 at t = 2 
with periodic boundary conditions and r/h = 0.8. 


Loo-ERROR 


L r ERROR 


N 

UN02 

TVD2 

'UN02 

TVD2 

20 

7.097 x 10- 3 

8.119 x 10" 2 

8.944 X 10“ 3 

6.778 x 10" 2 

40 

1.607 x 10 -3 

3.477 x 10~ 2 

2.044 x 10 -3 

2.033 x 10" 2 

80 

3.870 X 10" 4 

1.453 x 10“ 2 

4.926 x 10~ 4 

5.626 x 10~ 3 

160 

9.201 x 10“ 5 

5.975 x 10~ 3 

1.172 x 10- 4 

1.528 x 10~ 3 




























Table 3a. Numerical solutions of u, + uu x = 0 5 u(x,0) = sinirx, at t = 0.15 and r/h 
UN02 and TVD2 


1,0-ERROR 


L r ERROR 


N 

UN02 

TVD2 

UN02 

TVD2 

20 

1.890 x 10~ 2 

2.238 x 10" 2 

1.090 x 10" 2 

1.854 x 10' 2 

40 

5.712 x 10" 3 

1.054 x 10' 2 

3.034 x 10“ 3 

5.051 x 10~ 3 

80 

1.552 x 10“ 3 

4.422 x 10" 3 

7.771 x 10' 4 

1.340 x 10“ 3 

160 

3.985 x 10“ 4 

1.837 x 10' 3 

1.965 x 10“ 4 

3.621 x 10~ 4 


Table 3b. Same as Table 3a for FN02 and FVD2. 


, -ERROR 


L r ERROR 


N 

FN02 

FVD2 

JFN02 

FVD2 

20 

6.938 x 10“ 3 

2.091 x 10“ 2 

3.726 x 10' 3 

1.322 x 10' 2 

40 

1.959 x 10" 3 

1.054 x 10" 2 

8.869 x 10' 4 

3.835 x 10' 3 

80 

5.106 x 10~ 4 

4.424 x 10' 3 

2.163 x 10“ 4 

1.072 x 10~ 3 

160 

1.251 x 10~ 4 

1.837 x 10' 3 

5.270 x 10' 3 

2.946 x 10“ 4 






















Fisure la 


^i.curg 1 . Approximations of u = sin nx, -1 s x == 1, with N 
a) Piecewise-Parabolic interpolant Q(\,u). 
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PIECEWISE LINEAR TVD PIECEWISE LINEAR 



Figure lb Figure lc 

Figure 1 . Approximations of u = sinnx, -1 sx <1 with IJ = 10. 

(b) Pi ecewise-1 inear approximation (c) Piecewise-1 inear approximation L TVD2 

_ UN02 f \ 









Fi^re 3 a. for u given by (7.7) with 



If 








Figure Ha . Numerical solution of u t + u x = °> u ( x »°) Figure 4b , Same as Figure 4a using TVD2. 

defined by (7*7) at t = 2 with periodic 
boundary conditions N = 40 t / h = 0.8, 
using UN02. 








J S' 







Figurej6a. Solution to (7.8) with a = 2, 0=1, N = 20, Figure 6b. Same as Fig 6a using TVD2. 
CFL# = .5 , at t = 0.35 using UM02. ' — 






Figure 6C . Same as Fig. 6a using FN02. Figure 6D. Same as Fig. 6a using FVD2. 






Figure 7 « Numerical solution of density in a Riemann problem for Euler’s equations, 
(a) UN02 (b) TVD2 













’ Figure 9b . 

lfumerical solution of pressure in a Riemann problem for Eller's equations, (a) 0IW2. (b) T D2. 







